| Branch and Bound with Bernstein polynomials

> {-# LANGUAGE ViewPatterns, GeneralizedNewtypeDeriving,
>   MultiParamTypeClasses, FlexibleInstances #-}

> module OptBs (infIO) where

> import qualified BB
> import qualified Bs
> import qualified Data.List as List
> import qualified Aux
> import Aux (Hand(..))
> import qualified Iv
> import Iv (Iv)

> data (Fractional a, Ord a) => Prob a = Prob 
>   { probBs :: [Bs.Scaled a]
>   , probPath :: String
>   , probBox :: [Iv a] --TODO reconstruct from properly encoded probPath
>   , probBounds :: Iv a
>   } deriving Show

> prob :: (Ord a, Fractional a) => [Bs.Scaled a] -> String -> [Iv a] -> Prob a
> prob bs p x =
>   let l = foldl1 max (map Bs.lb bs) in
>   let bs' = filter (\b -> l <= Bs.ub b) bs in
>     Prob bs' p x $ Iv.iv l (Bs.lbubN (length x) bs')

> instance (Fractional a, Ord a) => BB.Box a (Prob a) where
>   dim = length . probBox

>   boxBounds = probBounds 

>   sharp p i = let s = map (Bs.sharp i) (probBs p) in
>               if all (Just L ==) s then Just L else
>               if all (Just R ==) s then Just R else Nothing

>   proj (Prob bs p x _) k h = prob (map (Bs.proj k h) bs) (p ++ show k ++ show h)
>                                     (Aux.remove k x)

>   pickDim (Prob bs p x _) = pickDir1 (length x) x bs

>   subdiv (Prob bs p x _) k =
>     let (x1,x2) = Iv.subdivBox k x in
>     (prob (map (\(Bs.Scaled b s) -> Bs.Scaled (Bs.subdivl k b) s) bs) (p ++ show k ++ "l") x1
>     ,prob (map (\(Bs.Scaled b s) -> Bs.Scaled (Bs.subdivr k b) s) bs) (p ++ show k ++ "r") x2)

>   boxInfo p _ = probPath p

> infIO :: (Fractional a, Ord a) => 
>      [Bs.Bs a] -- ^ A list of polynomials, whose maximum is to be minimized.
>   -> [Iv a]    -- ^ The domain.
>   -> a         -- ^ An \"epsilon\" that specifies how wide the returned interval may at most be.
>   -> Int       -- ^ Verbosity level
>   -> IO (Iv a, BB.Bx (Prob a)) -- ^ An interval enclosing the infimum and a minimizer.
> infIO bs box eps = BB.infIO (BB.box2bb $ prob (map Bs.mkScaled bs) "" box) eps

Strategies to select the dimension to be subdivided. These are
experimental.

> lowest :: Ord a => [(t, a)] -> (t, a)
> lowest = head . List.sortBy (\(_,l1) (_,l2) -> compare l1 l2)

> highest :: Ord a => [(t, a)] -> (t, a)
> highest = head . List.sortBy (\(_,l1) (_,l2) -> compare l2 l1)

> pickDir0 :: (Fractional a, Ord a) => Int -> [Iv a] -> [Bs.Scaled a] -> Int
> pickDir0 n x bs =
>   let (Bs.Scaled b _) = fst $ lowest $ Aux.map' (\b -> (b,Bs.lb b)) bs in
>     fst $ highest $
>     zipWith (\i x -> (i, let d = Bs.Scaled (Bs.diff i b) 1 in max (abs (Bs.lb d)) (abs (Bs.ub d)))) [0..(n-1)] x

> pickDir1 :: (Fractional a, Ord a) => Int -> [Iv a] -> [Bs.Scaled a] -> Int
> pickDir1 n x bs =
>   let (Bs.Scaled b _, _) = lowest $ map  (\b -> (b,Bs.lb b)) bs in
>     fst $ minBy (\(_,l1) (_,l2) -> l2 < l1) $
>     zipWith (\i x -> (i, let d = Bs.mkScaled $ Bs.diff i b in max (abs (Bs.lb d)) (abs (Bs.ub d)))) [0..(n-1)] x

> pickDir2 :: (Fractional a, Ord a) => Int -> [Iv a] -> [Bs.Scaled a] -> Int
> pickDir2 n x bs =
>   let (Bs.Scaled b _) = fst $ lowest $ map (\b -> (b,Bs.lb b)) bs in
>     fst $ highest $
>     zipWith (\i x -> (i, let d = Bs.mkScaled $ Bs.diff i b in max (abs (Bs.lb d)) (abs (Bs.ub d)) * Iv.width x)) [0..(n-1)] x

> pickDir3 :: (Fractional a, Ord a) => Int -> [Iv a] -> [Bs.Scaled a] -> Int
> pickDir3 n x bs =
>   fst $ highest 
>    (zip [0..(n-1)] x >>= \(i,xi) ->
>     bs >>= \(Bs.Scaled b _) ->
>       return (i, let d = Bs.mkScaled $ Bs.diff i b in max (abs $ Bs.lb d) (abs $ Bs.ub d) * Iv.width xi))

> minBy :: (a -> a -> Bool) -> [a] -> a
> minBy _ [] = error "minBy: empty list"
> minBy lt (x:xs) = minBy_ x xs
>   where minBy_ m [] = m
>         minBy_ m (x:xs) = minBy_ (if lt x m then x else m) xs
